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Abstract 

O . 

, In this paper, we propose a framework to solve a demand-supply optimization prob- 

^"*^ ' lem of long-term water resource allocation on a multi-connection reservoir network which, 

. , in two aspects, is different to the problem considered in previous works. First, while all 

' previous works consider a problem where each reservoir can transfer water to only one 

fixed reservoir, we consider a multi-connection network being constructed in Thailand in 
which each reservoir can transfer water to many reservoirs in one period of time. Second, 
, a demand-supply plan considered here is static, in contrast to a dynamic policy consid- 

(-^ ' ered in previous works. Moreover, in order to efficiently develop a long-term static plan, 

CZ3 , a severe loss (a risk) is taken into account, i.e. a risk occurs if the real amount of water 

i stored in each reservoir in each time period is less than what planned by the optimizer. 

The multi-connection function and the risk make the problem rather complex such that 
traditional stochastic dynamic programming and deterministic/heuristic approaches are 
inappropriate. Our framework is based on a novel convex programming formulation in 
which stochastic information can be naturally taken into account and an optimal solution 
jC' , is guaranteed to be found efficiently. Extensive experimental results show promising results 

■ of the framework. 
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1. Introduction 

X 

Water resource allocation in a water reservoir network is an important problem. In several 
decades, there has been intensive research on optimization of a reservoir policy in order 
to optimally control amount of water for each reservoir to match its demand in each time 
period (Murray & Yakowitz, 1979; Yakowitz, 1982; Archibald et al., 1997; Cai et al., 2001; 
Labadie, 2004; Cervellera et al., 2006; Castelletti et al., 2007; Wardlaw & Sharif, 1999; 
Tospornsampan et al., 2005; Reis et al., 2006; Jalali et al., 2007; Li & Wei, 2008; Wurbs, 
1993, 2005). In this paper, we focus on a different but related problem, namely, a demand- 
supply optimization problem of water resource allocation in a complex reservoir network. 

There are two major differences between ours and a problem of obtaining a dynamic 
policy considered in literatures. Firstly, existing literatures assume a simple network topol- 
ogy of reservoirs such that (1) a reservoir network is acyclic and (2) water released from a 
reservoir directly enters at most one other reservoir. This assumption is usually implicit in 
most of previous works but is stated explicitly by Archibald et al. (1997). However, cur- 
rently a more complex reservoir network has been being constructed in Petchburi province, 
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Figure 1: (Left) An example of the traditional reservoir topology assumed in existing works. 

(Right) An example of a multi-connection reservoir topology considered in this 
work. Solid lines denote possible transferring directions between two reservoirs. 
Dashed lines denote natural inflow and released water for consumers in each area. 



Thailand. In this multi- connection network, each reservoir can transfer water to many 
reservoirs in one period of time by using pipes which link between reservoirs. Figure 1 
illustrates the difference between the two topologies. Optimization in a multi-connection 
network is dimensionally much higher than a traditional network. 

Secondly, while the goal of existing works is to obtain a dynamic policy of water release, 
our goal here is to have a static long-term demand- supply plan. In a problem of obtaining a 
dynamic policy, it is assumed that demand is known for each reservoir in each time period. 
In our problem, demand (of each reservoir in each time period) itself will be optimized. 
Therefore, our problem can be considered as a prior step needed to be solved before solving 
a dynamic policy optimization problem. Moreover, in our problem, a risk management is 
also taken into account as there will be a risk (a severe loss) if natural inflow supply cannot 
meet the planned supply. More precisely, after achieving a long-term, e.g. year, demand- 
supply plan, the obtained plan will be declared to consumers such as farmers and industries. 
Consumers then integrate the declared supply to their year plans. In order to efficiently 
utilize their available water, they may invest some money, e.g. farmers may buy new cattle. 
Therefore, to prevent a futile loss of the investment, if real amount of inflow water is less 
than what is stated in the plan 1 , the government will have to pay an extra cost, here called 
a risk, of making extra water from somewhere outside the network, e.g. artificial rain, so 
that the amount of water will be as specified in the plan. 

As a problem is high-dimensional and not dynamic, it is inappropriate to apply stochas- 
tic dynamic programming (SDP) as in previous works (Cervellera et al., 2006; Castelletti 
et al, 2007; Murray & Yakowitz, 1979; Yakowitz, 1982; Archibald et al., 1997; Labadie, 
2004; Wurbs, 2005). Beside SDP, deterministic and heuristic approaches such as linear pro- 
gramming, quadratic programming, genetic algorithm and simulated annealing have also 
been previously employed (Cai et al., 2001; Wardlaw & Sharif, 1999; Tospornsampan et al., 
2005; Li & Wei, 2008; Jalali et al., 2007; Labadie, 2004; Wurbs, 2005). Nevertheless, these 
approaches do not directly support the use of stochastic information which is fundamental 
to reservoir network optimization. In these works, the stochastic information is usually 



1. Due to the stochastic nature of the rainfall and run-off phenomena. 
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implicitly employed. For examples, given a probability distribution of inflow for each reser- 
voir in each time period, the expected value of inflow can be calculated and input to a 
deterministic (or heuristic) method. 

In this paper, we propose to apply the framework of convex programming to solve the 
problem. The nowadays technology of convex programming is very efficient like least-square 
and linear programming, and, with some efforts, many non-linear problems can be reduced 
to convex programming formulations (Boyd & Vandenberghe, 2004). Similar to linear 
programming, although convex programming is efficient and appropriate to solve the static 
problem, a direct application of convex programming does not naturally support the use of 
stochastic information. Here, we present a novel convex programming formulation which is 
able to directly take the stochastic information into account and yield a reasonable output. 

The remaining of the paper is structured as follows. In Section 2, a mathematical for- 
mulation of our problem will be given. In Section 3, related frameworks and our framework 
which is proposed to solved the problem will be explained. After that, experimental results 
and discussions are provided subsequently. 



2. Problem Specification 

Our problem specification here is similar to that of Cervellera et al. (2006). However, 
Cervellera et al. considered only a simple network in Figure 1 (left) and did not consider a 
static plan with risk management. Let TV" be the number of reservoirs in the network and 
T be the number of periods considered in the plan, e.g. one year. Let indices n, m run 
through 1, ...,N and t through 1,...,T. At a period of t a reservoir n can transfer some 
water to another reservoir m with amount 

< q l nm < Qnm, (1) 

where Q nm is the maximum pumpage capacity in one period between the two reservoirs. 
There is also a pumpage cost 

where C nm (-) is a convex and non-decreasing function. Aside from inputs from other reser- 
voirs, each reservoir will also have its natural inflow input, 

xi ~ Pr(xi), (3) 

which is a random variable with respect to a known probability distribution Pr(x n ). For 
each month, there is a profit from releasing water to consumers 

<£(</»), (4) 

where G n (-) is a concave and non-decreasing function, and g n is an amount of released water 
to consumers. Let v n be a volume (after transferring q nm and releasing g^) of each reservoir 
on each month. There is a maximum volume M n for each reservoir. If v n < M n , according 
to the traditional rule of conservation, v n will depend on g n , x n and q nm - These relations 
can be concluded by the following state equation of a volume: 



vi = min < 



M n, vi 1 - g l n + X l n + Qmn ~ Yl ^ \ • ( 5 ) 
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And the amount of released and transferred water cannot exceed the existing volume; hence, 

0< 5 *+^^ m <<- 1 . (6) 

n 

As noted by Cervellera et al. (2006), Eq. (6) is conservative since we do not include both 
(uncertain) natural and transferred inflow at period t. In the problem, the initial volume 
of each reservoir v® has to be specified. Moreover, there are volume constraints at the end 
of the optimization period: 

Vl > V n . (7) 

Notice that, unlike previous works, there are two notions of sending water from a reservoir, 
i.e. water transfer q i nm and water release g l n . In contrast to water release which is for 
consumption, water transfer between two reservoirs requires some amount of energy for a 
pumpage. Hence, there is a cost for water transfer, but a profit for water release. 

Finally, we take a risk into account. Since the purpose of a long-term plan considered 
in this paper is to optimize the demand-supply, after obtaining g l n from an optimizer, the 
resulted g l n will be declared to consumers such as farmers and industries. Hence, g l n will be 
their target demand (or, a target supply for the government) since consumers will integrate 
g l n to their long-term plans. In order to efficiently utilize their available water, consumers 
may invest some money, e.g. farmers may buy new cattle. As a result of a stochastic nature 
of x l n stated in Eq.(3), realizable amount of release water g l n may different to its target g l n . 
If g l n < g l n , there will be a severe extra cost, or a risk, of acquiring extra water, e.g. by 
using artificial rain, to satisfy the target demand for consumers: 

RWn-tfn) (8) 

where i?^(-) is a convex, non-decreasing function and is zero if g l n > g l n . 

The goal of this optimization problem is to maximize a profit (or minimize a loss) subject 
to all constraints described above where optimized variables are {^m,^}- 

3. The Proposed Framework 

In this section, we provide our convex optimization framework to solve the problem. How- 
ever, at first, we motivate the necessary of our framework by illustrating the difficulties of 
using traditional methods to solve the problem. 

3.1 The Need of a New Problem Formulation 

From the previous section, it can be observed that our optimization problem is similar to 
traditional problem formulations of dynamic policy optimization; one may temp to think 
that methods proposed in previous works can be employed to solve our problem as well. 
The purpose of this subsection is to explain reasons that existing methods are in fact 
inappropriate for our current problem, and thus a new method is needed. 

We begin by considering the most-popular stochastic dynamic programming (SDP) 
method (Murray & Yakowitz, 1979; Yakowitz, 1982; Archibald et al, 1997; Labadie, 2004; 
Wurbs, 2005; Cervellera et al., 2006; Castelletti et al., 2007). There are some difficulties to 
solve the problem by using SDP due to its nature of dynamic closed-loop control. Firstly, 
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Figure 2: The traditional SDP formulation of the multi-connection reservoir optimization 
problem. This formulation cannot take a risk into account and is intractable. 



since a long-term static plan is our target, all {^ m , gl l } t _ 1 must be obtained at the be- 
ginning period t = 1. However, for any t > 2, a policy obtained from SDP will provide 
us optimal values of <2^ m , 5^ only when the state variables {t^ -1 , v^}, which are not known 
at t = 1, are provided (Russell & Norvig, 2003, Chapter 17). Secondly, the risk defined 
in Eq. (8) cannot be straightforwardly taken into account since, by its definition, a risk 
occurs because a decision has to be made before knowing the real situation (i.e. a decision 
is made before knowing the future values of state variables); nevertheless, a policy obtained 
from SDP would provide a decision only after the real situation is known (i.e. a decision 
will be made in the future). These difficulties may be crudely solved by first implementing 
traditional closed-loop SDP (without the risk) as illustrated in Figure 2, and then employ 
the Viterbi algorithm (Russell & Norvig, 2003, Chapter 15) to find the most-likely sequence 
of {q nm ,9n\- Therefore, a series of dynamic controls can be converted to a static plan by 
this hybrid SDP- Viterbi algorithm. To take the risk into account, some algorithms may 
be used to further hybridize the SDP- Viterbi algorithm. However, this complicated hybrid 
method of converting sequential dynamic-controls to a static plan with risk management is 
beyond the scope of this paper. 

In fact, we note that there are a third difficulty arising when using SDP since the 
SDP formulation illustrated Figure 2 is not only stochastic and non-linear, but also high- 
dimensional. The number of variables grows quadratically with respect to N (due to the 
q nm variables); thus, dimensions of the problem is relatively much higher than problems 
considered in previous works. Due to the curse of dimensionality (Bellman, 1954; Larson, 
1968; Bertsekas & Tsitsiklis, 1996) of discretizing the variables, this high-dimensionality 
actually prevents the use of SDP in general; even the most efficient SDP approximation, 
the neuro- dynamic programming method (Bertsekas & Tsitsiklis, 1996), known in literatures 
can cope with problems of only 30 dimensions (Cervellera et al., 2006). 

Beside SDP, some existing works proposed to apply deterministic and heuristic ap- 
proaches such as linear programming (LP), quadratic programming (QP), genetic algorithm 
and simulated annealing for conventional reservoir networks (Cai et al., 2001; Wardlaw &i 
Sharif, 1999; Tospornsampan et al., 2005; Li & Wei, 2008; Jalali et al., 2007; Labadie, 2004; 
Wurbs, 1993, 2005). Nevertheless, the most important limitation of these methods is that 
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they cannot directly take the stochastic information and the risk into account 2 . Usually, 
uncertainty of inflow is ignored, and x l n is simply replaced by a deterministic value such as 
its mean x\. Henceforth, a deterministic approach with these simple deterministic value 
inputs will be called a traditional (or conventional) deterministic method. Since a heuristic 
approach does not guarantee an optimal solution, the traditional deterministic approach is 
the most suitable method among existing approaches for solving the problem considered 
here. 

Finally, we note that in the literatures of convex optimization itself, there is a method- 
ology called robust programming (Boyd &; Vandenberghe, 2004, Chapter 4) which indeed 
can efficiently utilize the stochastic information. Nonetheless, this robust programming is in 
fact a worst-case optimization method where the worst situation will be solved. Specifically, 
in our reservoir problem, robust programming will consider the situation where minimum 
inflow occurs for each reservoir in each time period. This worst case situation indeed 
rarely occurs, and hence robust programming is also inappropriate for our problem. In the 
next section, we present a novel convex optimization formulation which can corporate the 
stochastic information in a more-efficient way. 

3.2 Convex Programming Formulation 

Here, we give a new optimization formulation to appropriately solve the problem stated in 
Section 2 using the convex programming framework (Boyd & Vandenberghe, 2004). Note 
that in contrast to stochastic dynamic programs, high-dimensional convex programs which 
are deterministic can be solved efficiently. Moreover, the long-term plan where {<7^ m , 
is needed at the beginning period can be easily achieved. 

The most important step here is to take the stochastic inflow information into this 
deterministic framework. Our main idea is, instead of treating stochastic inflow variables 
as random variables or constants like previous works, the stochastic variables x l n themselves 
will be optimized, as our "best predictions" of inflows. Then, in order to take the risk into 
account, we propose to consider the average risk of Eq.(8) over all possible actual inflow. 
The best inflow predictions will be optimized so that they will balance between profits 
and risks. For instances, if the risk function R(-) is severe, the best predictions will be 
conservative such that they will smaller than the mean, i.e. x n < x^; this small- value 
prediction of x n will limit a predicted volume v n , and hence also limit amount of water 
transfer q l nm: i.e., if the risk is costly and we are unsure of being have abundant of water, 
we should not transfer the water to other reservoir carelessly; in contrast, if R(-) return low 
values, the best predictions will be aggressive such that they will be smaller than the mean, 
i.e. x n > x n . 

To formally explain our framework, we first note that if we are able to guarantee that 
for each reservoir and each time period v l n < M n , then the state equation simplifies to 

vl = v'n 1 -g'n + xi+Yl q ™ ~ q nm> ■ ( 9 ) 



2. Because the risk is defined based on uncertainty, if uncertainty itself cannot be taken into account, the 
risk also cannot be incorporated. 
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Figure 3: Summary of the proposed optimization problem. 



It will be explained below that our problem formulation guarantees to satisfy the condition 
v\ < M n . Now, suppose an actual (future observed) inflow is x\. Define g^ as the realizable 
amount of water released by a reservoir n at a period t (see Eq.(8)) 

g t n =9 t n + V t rT 1 -Vn 1 , (10) 
where v n is as an actual volume defined recursively as 

and 4 = u„. Eq.(lO) states that the realizable amount of release from the reservoir n at 
period t is the target release g l n plus the difference between the actual and target volumes 

4" 1 - «n -1 - % E q-( 10 )> we have that 9n ~9n = 4" 1 and 4" 1 -9n = 4 _1 Using 

these relations and subtract Eq.(ll) with Eq.(9), then, the expected risk E^ \j& n (4 — 4)] 
can be stated as 

E x [J?* (& - 0* )] = Ex [R l n (x 1 - 1 - x 1 - 1 )] ■ (12) 

For any t, to (approximately) calculate Ej fa^ — 4)1 j we discretize the domain of 
Pr(x t n ) into a finite domain set, X\: 

E £ (4 - 4)] « ^ Pr(4)i?* (4 - 4) . (13) 

Note that in practice the domain of Pr{x t n ) is originally discrete since we are able to collect 
only finite statistics of x\ and thus Eq.(13) can be actually exact. Note further that since 
Rn(-) is convex, &x Pr(x t n )R t n (x^ — x^) is a l so convex. 

Our optimization formulation is summarized in Figure 3. This maximization problem 
formulation consists of a concave objective function and linear constraints, and thus the 
problem can be efficiently solved by the convex programming framework (Boyd & Vanden- 
berghe, 2004). The initial values of v® has to be provided. 

Note that in order to make our problem formulation convex and efficient, some con- 
straints in Section 2 are changed and new constraints are added. First, the constraint 
"min(Af^) < x f n < max(Af^)" is added so that x l n will take only realizable values. Next, 
the constraint of Eq.(5) is simplified by introducing a new cost term "i* 1 * max(^ — M n , 0)". 



7 



Chatpatanasiri & Sriburi 



Table 1: Summary of differences between the three frameworks. 





Traditional SDP 


Traditional LP/QP 


Our Convex Program 


Type 


Closed-Loop Control 


Open-Loop Optimization 


Open-Loop Optimization 


Policy 


Dynamic: Real-Time Decision 


Static: Long-Term Plan 


Static: Long-Term Plan 


Running Time 


Intractable 


Tractable 


Tractable 


Stochastic 


Included 


Not Included 


Included 


Objective 


Non-linear 


Linear/Quadratic 


Concave (Maximization) 


Risk 


Not Included 


Not Included 


Included 


Solution 


Sub-Optimal 


Optimal 


Optimal 


X* 


Random Variables 


Constants 


Optimized Variables 


gl 


Control Variables 


Optimized Variables 


Optimized Variables 




Control Variables 


Optimized Variables 


Optimized Variables 



Here, is a big constant, and by setting large enough, this will guarantee that v n < M n . 
Note that to obtain a standard objective function form, this term can be easily replaced by 
a slack variable. 

Table 1 summarizes the differences between our framework and the other frameworks. In 
the table, we do not consider heuristic methods such as genetic algorithm or simulated an- 
nealing. The features of heuristic methods are similar to those of deterministic methods: the 
major difference is that heuristic methods can handle non- linear objective functions. Never- 
theless, an ability to handle a non-linear function comes up with a certain tradeoff: heuristic 
approaches are either sub-optimal or intractable. Also heuristic methods usually were tested 
with problems with convex objective functions (Wardlaw & Sharif, 1999; Tospornsampan 
et al., 2005; Jalali et al., 2007; Li & Wei, 2008). Therefore, our framework presented here 
can substitute these works with the main advantage that our framework guarantees to find 
an optimal solution with a polynomial running time. We note also that our framework in 
fact generalizes traditional uses of LP and QP: if we specify the Pr{x t n = x') = 1 for some x' 
(such as x^), then by the constraint min(A?*) <= x l n <= max(A"*), the best prediction will 
be x', and the risk term R{x l n — x^) will become a constant and can be removed from the 
optimization problem; for example, in cases that G(-) and C(-) are linearizable, the resulted 
program is simply an instance of traditional LP applied to a multi-connection reservoir 
network, see e.g. Wurbs (2005). 

4. Experiments 

As demonstrated in the previous section (see also Table 1), a deterministic method is the 
most appropriate candidate among existing methods which can be implemented to solve 
our problem; hence, in this section, we evaluate our method and compare it to a standard 
deterministic method, matlab with the optimization interfaces from Yalmip (Loefberg, 
2004) is used for all implementations. For each simulated situation, the following process is 
employed: 

1. All functions, C* (•), G^(-) and i?^(-)> all constants, Q nm and M n , and all probability 
distributions Pr(x^) are input to our method (or a deterministic method). The long-term 
demand-supply plan of {<^ m , g^} is then obtained for each method. 
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Figure 4: (Left) A simple 2-reservoir network. (Right) The Ang-Puang multi-connection 
reservoir network being constructed in Petchburi province, Thailand. The maxi- 
mum volume of each reservoir is indicated. Subscripts denote orders of reservoirs. 




2 4 6 -3 -2 -1 1 2 3 D 

Figure 5: The C^(-), R l n {-) and G^(-) employed in the Section 4.1 and Section 4.2. 

2. A sequence of inflow is generated {x l n } where each x l n ~ Pr{x t n ). is then calcu- 
lated according to Eq.(10) and Eq.(ll). 

3. The total profit of each method will calculated from: 

(fl'n) ~~ X^Cnm (q^m) — ^-^n (<?n ~~ 9n) ■ 



4. Steps 2. and 3. are repeated 100 times to calculate the average performance of each 
algorithm 3 . 



4.1 Simple Network 

To intuitively understand results obtained from each algorithm, it is best to start from a a 
simple reservoir network as illustrated in Figure 4 (Left). Here, C^(-), -R^(-) and G^(-) are 
illustrated in Figure 5 where several values of their slopes and £ D' (shown in the figure of 
Gr^(-)) will be tested. L D' can be interpreted as the "maximum demand" for each reservoir 
at a certain period. In this experiment, these functions are the same for all n and t. Since 
these functions are linearizable, a deterministic method considered in this example is LP. 
The following values are specified: T = 3, N = 2, Mi = M 2 = 10, Q12 = Q21 = 5, and 

3. Note that Step 1. is deterministic and not depend on a real sequence of inflow {x^}- 
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Figure 6: [The Top 4 Figures]: (Up-Left) The probability distributions of inflow used in 
the simple network. With the exception of Pr(x\) and Pr^), the distributions 
are the same for all n and t. (Up-Right) Performance summary: "(Average) 
Total Profit = Release Profit - Transfer Cost - (Average) Risk" with standard 
deviations of "Total Profit" and "Risk". The C,R,G slopes and the D value of 
Figure 5 are also shown on the top. (Down) {x^} and {fl^} are shown for each 
reservoir. [The Bottom 4 Figures]: Another case of the the simple network. 



< 2 



V\ = V% = 1. Figure 6 illustrates two cases of this simple network and the 
corresponding results obtained from traditional LP and our method. The traditional LP 
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Figure 7: The geography of the Ang-Puang multi-connection reservoir network. "In 
Progress" denotes a pipe which is being constructed or is planned to be con- 
structed. 



implementation can be achieved by simply replacing each uncertain x\ with a deterministic 
value; x l n is used here. 

The two cases here illustrate the same simple situation which can be intuitively explained 
as follows. The first reservoir has a good probability to have more than usual inflow at t = 1, 
and the second reservoir will surely do not have any inflow at t = 2. Therefore, clearly, 
the first reservoir should send water to the second reservoir on the second period, but the 
question is: how much should be transferred? The results obtained from traditional and 
our methods can be used as answers to this question. From the simulations of these two 
cases (and other several cases not shown here due to the space limitation), our method 
always significantly outperforms traditional LP. The main reason is because traditional LP 
inefficiently takes stochastic information into account. 

The top 4 figures of Figure 6 illustrate a case where the risk cost is high, and our method 
becomes more conservative about its predictions of inflow. Even though our method cannot 
make a high average released profit as that of traditional LP, our method does provide a 
much lower average risk cost, and therefore results in a higher average total profit (5.49) 
than that of traditional LP (4.91). The bottom 4 figures of Figure 6 illustrate a second case 
where a risk plus a transferred cost is not severe and, with high probability, x\ will be very 
high. In this case, our method results in aggressive inflow prediction of x\ which in turn 
leads to a higher release profit and a higher average total profit (7.99) compared to that of 
traditional LP (7.88). 
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Figure 8: (Up-Left) Inflow probability distributions of the three big reservoirs for t = 1,2. 

(Up- Right) Inflow probability distributions of the other reservoirs for t = 1,2. 
(Down-Left) vq and V n for each reservoir. (Down-Right) Summary of the profit 
for each method. 



4.2 Complex Network 

Currently, by royal initiating of his Majesty the King of Thailand, a multi-connection 
reservoir network, namely Ang-Puang, has been being constructed in Petchburi province, 
Thailand 4 . The geography of the network is shown in Figure 7, and the network topol- 
ogy, together with the capacity of each reservoir (in metric tonne), is shown in Figure 4 
(Right). For simplicity, we will mention each reservoir by its order number, e.g. reservoir-1, 
reservoir-2, etc., instead of its name. Since the network is currently being constructed and 
all the historical data are being collected, in this experiment we are able to employ only 
artificial data. The artificial data employed here is in the same spirit as the standard 4 and 
10 reservoir problems for a traditional network (Larson, 1968; Murray &i Yakowitz, 1979; 
Wardlaw & Sharif, 1999; Tospornsampan et al., 2005; Jalali et al., 2007). Nevertheless, our 
artificial data contains also stochastic information while there is none in the 4 and 10 reser- 
voir problems. Here, C*(-), R l n {-) and G l n {-) are the same as the previous case, illustrated in 
Figure 5. The following values are specified: T = 6, N = 8, Q nm = 2.5, G n = 1, C t nm = 0.25 
and R n = 2.5 for all n, m, t. Other parameters are shown in Figure 8 and Figure 9, i.e. 
Figure 8 illustrates Pr(x t n ) for t = 1,2 and v® = v® = V\ = V2, and Figure 9 shows the 
values of "maximum demand" D specified in G n {-) for each reservoir in each period. 

This problem setting is a simplified version of a real situation normally occured in 
Thailand where t is a two- month period and 1,...,T is thus one year. A long-term plan 
is determined just before the rainy season; thus, only for t = 1,2 where x l n will have 
probabilities to be significantly more than zero, and almost zero in the t = 3, 6 periods. 
To save space, Pr(x t n ) for t = 3, ...,6 (which their realizable values are around zeros) are 

4. http:/ /www. huaysaicenter.org/water.php 
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Figure 9: Inflow prediction x l n , released volume g l n , and D^, for each method and for each 
reservoir. 

not shown. The intuition of this problem setting is that, at the first two periods t = 1,2, 
because of the huge amount of rain, all reservoirs will, with high probability, nearly attain 
their maximum capacities. However, only the big reservoirs n = 1,4,8 will have enough 
water for later periods and are responsible to transfer water to the other small reservoirs. 

The total profit of the two methods are shown in Figure 8 (Down-Right). Note that, 
for both methods, standard deviations are almost zeros, and thus our method significantly 
provides more profit than the conventional deterministic method. As can be seen from 
Figure 9, in this situation, our method usually becomes more conservative about amount 
of inflow than the conventional method. As a result, our method cannot guarantee water 
supply in reservoir-6 and reservoir-7; in contrast, the plan of conventional LP decide to 
provide almost enough water for every reservoir in every period. Nevertheless, there are 
often cases that the real inflow x* for each reservoir is not as high as expected from x l n 
of the plan of conventional LP; thus, a high risk cost has often to be paid for extra water 
requirement. 

It is also important to note that, in the simulations, the time-complexity of our method 
increases only slightly from traditional LP. By using MATLAB's "linprog" function with the 



13 



Chatpatanasiri & Sriburi 




70 r 



60 - 



Sensitivity Analysis of "G" 



* - Tradiitional LP 
-o- Our Method 



1.2 1.4 1.6 1.8 2 

G 



35 r 



Sensitivity Analysis of "vO" 



30 - ^x-x-x* * 
1 25 - 

Q. 

o 20 - 

15 - 



- * - Tradiitional LP 
-o- Our Method 



\ e 



10^ 
o 



0.2 0.4 



0.6 0.8 
vO (% of M) 



Figure 10: Sensitivity analysis of the two methods. 



"largescale" option. Our method converges in 18 iterations (4.52 seconds) where traditional 
LP converges in 17 iterations (4.28 seconds). 

4.3 Sensitivity Analysis 

In this section, we demonstrate the sensitivity of each method with respect to four param- 
eters of the complex- network problem explained in Section 4.2. Figure 10 shows results 
of the sensitivity analyzes. In each analysis, all the parameters are the same as those of 
the experiment in Section 4.2, except the analyzed parameter written in the title of each 
figure. Each figure shows the average total profit of each method in each parameter setting. 
The standard deviation of the average is also drawn; nevertheless, the standard deviation 
of each point is very small and is difficult to observe in every figure. This indicates that 
the statistical confidences of the results shown here are very high. Note that our method 
always significantly outperforms the existing method in all cases. 

Figure 10 (Up-Left) illustrates the case where the parameter 'C^ m ' is varied. We can 
see that in the beginning where C t nm < 1, the total profits of the two methods decrease 
about linearly with respect to the increase of C i nm . After C l nm > 1, the total profits do not 
change further since the transfer cost becomes more than the release profit G^ m , and thus 
no water transfer occurs after ^> 1. 

Figure 10 (Up-Right) illustrates the case where the parameter l R t n ' is varied. We can see 
that the average total profit of traditional LP decreases about linearly with respect to the 
increase of R l n . In contrast, the average total profit of our method decreases in a slower rate. 
This is because our method is able to efficiently take the risk and the stochastic information 
into account while traditional LP cannot. As a result, our method is less sensitive to the 
change of than traditional LP. 

Figure 10 (Down-Left) illustrates the case where the parameter 'G^' is varied. In this 
case, the average total profits of the two methods increase about linearly with respect to 
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G n . In this case, like other cases, our method provides significantly higher total profit than 
traditional LP although it is difficult to see in the figure due to the large scale of the y-axis. 

Figure 10 (Down-Right) illustrates the case where the parameter i v^ l = V n ' is varied for 
each reservoir. In this case, the is varied from 5% of M n to 100% of M n . As a result, the 
sensitivity of each method is about the same. Note that as the value of v ® = V n becomes 
near M n , the total profits decrease since the constraint >V n forces the optimizer not to 
use any water at later periods. 

5. Discussions 

As we explained earlier, a deterministic method seems to be the most appropriate among 
previously existing methods for the problem of demand-supply optimization in a multi- 
connection network considered in this paper. From the experiments, we have shown that our 
new optimization framework gives promising results by always outperforming the existing 
deterministic method in all experiment settings. 

There are some issues which should be considered in practical uses of our method. 
Firstly, remember that our method requires to discretize Pr(x^). This discretization is not 
an important limitation of our method since, in practice, all statistics from historical data 
are in fact originally discrete. Hence, the original statistical distributions can be applied to 
our method directly. Nevertheless, an application of a smoothing method might be useful 
to obtain simplified probability distributions which can reduce an optimization time of our 
method. 

Secondly, we note that although the class of convex programming is efficiently solvable, 
the class of linear programming is still the most preference in term of time efficiency In fact, 
theoretically, a convex function can be linearized which an arbitrary degree of approximation 
accuracy. Therefore, in our opinion, one potential future research direction is to apply an 
efficient linearization method from the computer vision community to linearize a convex 
objective function (Salotti, 2001; Kolesnikov & Franti, 2003), and find the best tradeoff 
between the time efficiency and the approximation accuracy. 

Finally, it is important to note that although our method provides a static long-term 
plan, our method can be used to dynamically update a plan as well (by re-running the 
method in the beginning of each time period), provided that the length of t is not shorter 
than a running time of our method. In all experiments mentioned in Section 4, the opti- 
mization process of our method usually accomplished within 1 minute; therefore, for any 
long-period plans, e.g. t is a week or a month, our method can also provide a dynamic plan. 
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